!   _______________________________________________________________________________
!    _______/\\\\\\\\\_______/\\\\\\\\\______/\\\\\\\\\\\\\\__/\\\\\\\\\\\\\\\______
!     _____/\\\////////_____/\\\////////_____\/\\\//////////__\/\/////////////________
!      ___/\\\/____________/\\\/______________\/\\\____________\/\\\________/\\\\_______
!       __/\\\____________/\\\_________________\/\\\\\\\\\\\\\__\/\\\_______\////________
!        _\/\\\___________\/\\\_________________\/\\\/////////___\/\\\________/\\\\_______
!         _\//\\\__________\//\\\________________\/\\\____________\/\\\_______\////________
!          __\///\\\_________\///\\\______________\/\\\____________\/\\\______/\\\\_________
!           ____\////\\\\\\\\\____\////\\\\\\\\\___\/\\\____________\/\\\\\\\\\////__________
!            _______\/////////_______\/////////_____\///_____________\/////////_______________
!            __________________________________________________________________________________
!----------------------------------------------------------------------------------------------
!>  subroutine input_files_cgns
!>  last edit 2016-05-24
!>  last edit by liuxz
!----------------------------------------------------------------------------------------------
	subroutine sta_aero(hight,re,ma,cref)

        real*8 :: temp,h,cre,re,ma
        real*8 :: t,p,rho,a,g,miu

        real*8 :: t_sl=2.8815e2,p_sl=1.01325e5,rho_sl=1.2250
        real*8 :: g_sl=9.80665, a_sl=3.40294

        real*8 :: r0=6.356766e6

        h=hight/(1+hight/r0)/1000

        if(0<=hight.and.hight<=11.0191e3) then

            temp=1-h/44.3308
            t=288.15*temp
            p=temp**(5.2559)*p_sl
            rho=temp**(4.2559)*rho_sl

        else if(11.0191e3<hight.and.hight<=20.0631e3) then

            temp=exp((14.9647-h)/6.3416)
            t=216.65
            p=0.11953*temp*p_sl
            rho=0.15898*temp*rho_sl

        else if(20.0631e3<hight.and.hight<=32.1619e3) then

            temp=1+(h-24.9021)/221.552
            t=221.552*temp
            p=2.5158e-2*temp**(-34.1629)*p_sl
            rho=3.2722e-2*temp**(-35.1629)*rho_sl

        else if(32.1619e3<hight.and.hight<=47.3501e3) then

            temp=1+(h-39.7499)/89.4104
            t=250.350*temp
            p=2.8338e-3*temp**(-12.2011)*p_sl
            rho=3.2618e-3*temp**(-13.2011)*rho_sl

        else if(47.3501e3<hight.and.hight<=51.4125e3) then

            temp=exp((48.6252-h)/7.9223)
            t=270.650
            p=8.9155e-4*temp*p_sl
            rho=9.4920e-4*temp*rho_sl

        else if(51.4125e3<hight.and.hight<=71.8020e3) then

            temp=1-(h-59.4390)/88.2218
            t=247.021*temp
            p=2.1671e-4*temp**(12.2011)*p_sl
            rho=2.528e-4*temp**(11.2011)*rho_sl

        else if(71.8020e3<hight.and.hight<86e3) then

            temp=1-(h-78.0303)/100.2950
            t=200.590*temp
            p=1.2274e-5*temp**(17.0816)*p_sl
            rho=1.7632e-5*temp**(16.0816)*rho_sl

        else if(86e3<=hight.and.hight<=91e3) then

            temp=exp((87.2848-h)/5.47)
            t=186.870
            p=(2.2730+1.042e-3*h)*10e-6*temp*p_sl
            rho=3.6411e-6*temp*rho_sl

        end if

        a=sqrt(1.4*287*t)
        g=9.80665/(1+h/6.356766e3)**2
        miu=(1.458e-6*t**1.5)/(t+110.4)
        v=ma*a
        re=(rho*v*cref)/miu

		return

	end subroutine sta_aero
